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CO ■ Abstract. - We present a simple mesoscale field theory inspired by rate-independent plasticity 

that reflects the symmetry of the deformation process. We parameterize the plastic deformation 
by a scalar field which evolves with loading. The evolution equation for that field has the 
form of a Hamilton-Jacobi equation which gives rise to cusp- singularity formation. These 
cusps introduce irreversibilities analogous to those seen in plastic deformation of real materials: 
we observe a yield stress, work hardening, reversibility under unloading, and cell boundary 
formation. 
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We call it plasticity when materials yield irreversibly at large external stresses. Macroscop- 
ically, plasticity is associated with three qualitative phenomena. To a good approximation, 
there is a threshold called the yield stress below which the deformation is reversible (see 
Fig. [l]). A material pushed beyond its yield stress exhibits work hardening, through which 
the yield stress increases to match the maximum applied stress. Finally, the deformed crystal 
develops patterns, such as the cell structures observed in foe metals ^,|2[^. While much is 
known about all three phenomena, a quantitative understanding based on mesoscopic con- 
tinuum theory would be welcome — especially if it connects to microscopic properties of the 
atomic interactions in the material. In this paper, we will discuss a simple scalar field theory 
which naturally exhibits these three key features of plasticity. We do not claim to model in 
details plasticity in real materials but we believe that a continuum description of plasticity 
should share the key feature of our model equation: The transition between reversible and 
irreversible deformation is generated by singularities which occur at finite stress. 
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Fig. 1 - Generic stress-strain relation. For stress a below the critical yield stress a yo the deformation 
is elastic and reversible. Further deformation results in plastic deformation and work hardening, 
i.e. the yield point (bullet) moves to higher stress. The strain e is partly due to reversible elastic 
deformation e c and partly due to irreversible plastic deformation e p . 



There are a vast number of approaches to modeling plastic flow in metals — quantum, 
atomistic, motion of single dislocations, motion of many dislocations, continuum theories of 
dislocation densities, slip system theories, work hardening theories .... Reviewing even only 
the basic ideas of these models is beyond the scope of this paper. Three of them, however, 
directly inspired our approach. (1) Theories based on the differential geometry of the Burgers 
vector density torsion tensor These elegant mathematical descriptions of the state of 

the material need to be supplemented by a similarly sophisticated dynamical evolution law — 
especially for the non-equilibrium problem of plastic flow. The torsion tensor theories are 
typically dismissed by the engineering community because they ignore the large majority of 
geometrically unnecessary dislocations which cancel out in the macroscopic Burgers vector 
density. (2) Macroscale engineering theories of plasticity, and in particular the recent strain 
gradient theories 0. Their use of symmetry to constrain the form of the evolution laws is 
echoed in our approach. Their explicit incorporation of a yield surface is clearly appropriate on 
the macroscale, but is unsatisfactory for a condensed-matter physicist. In our model, the yield 
surface emerges naturally from the evolution laws through the development of cell structures 
on the mesoscale. (3) Theories which attempt to coarse-grain the complex rearrangement 
dynamics of atoms or dislocations || to develop continuum plasticity theories. Any such 
description will introduce a field which describes the local state of the material as it evolves 
during deformation. Even though the microscopic dynamics of dislocations is understood to 
a large extent, coarse-graining the microscopic dynamics has not been successful up to now. 

We are going to study the following evolution equation for a scalar field which parameterizes 
rate-independent plastic deformation 

We use Einstein's convention, summing over repeated indices. The plastic deformation is 
assumed not to depend on volume changes and the evolution depends only on the deviatoric 
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stress Sij = cr^ — | er^ Sij, with er^ the full stress tensor. There are two motivations to study 
this equation in the context of plasticity theory. 

First, we can derive this evolution equation as one of the leading terms in a long- wavelength 
expansion of a general local evolution equation for ^ in which respects the symmetries and 
the rate invariance of the loading process. However, there is a priori a term of the same order 
ViVj'J' which turns out to be highly singular since it it leads to finite time divergences 
of the field Although rate independent (creep) terms such as V 2 \l/ (and higher order 
derivatives) should regularize these divergencies, we prefer to discard these singular terms 
altogether in the present discussion and focus on eq. (^) which, as we now show, already 
contains an rich phenomenology. We will indeed show later that eq. (Q) can be transformed 
into a Hamilton- Jacobi equation which is known to form finite time singularities. These 
singularities are the second reason to study this equation since they can be related to the onset 
of irreversibility and to dislocation structure formation. Second, Hamilton- Jacobi equations 
are also related to conservation laws as studied in the context of traffic-jams and at least on 
a naive level there is a similarity between traffic-jams and dislocation entanglement. 

We will focus on so-called proportional loading paths (such as those used in shear and 
tension tests, where the loading direction does not change) where the stress Sij(t) — a(t) Sij. 
Thus in eq. ([!]) we change variables from time t to stress amplitude a 

<9\I> 

— =S y V,*V i *. (2) 

In the form of eq. (||), we see the key challenge in formulating laws of rate-independent 
plasticity: the equations appear to be manifestly reversible. Increasing and then reducing a 
will naively leave the material in the original state. The engineers bypass this problem by 
formulating their theories not in terms of order parameter fields, but directly in terms of a 
yield surface (corresponding to step functions 0(er — a y ) in the equations of motion). However, 
eq. (^|) is a Hamilton- Jacobi equation (closely related to the multidimensional anisotropic 
Burgers equation, see JlOj] ) which develops cusp-shaped singularities at finite a even for smooth 
initial conditions. At these singularities, as numerically illustrated below, information and thus 
reversibility is lost. The stress at which the first singularities form has many features of a 
yield surface as we will show later. 

After the formation of singularities the solution of eq. (^) is not unique and one has to find 
a weak solution which regularizes the singularities. As mentioned above, the assumption of 
a rate-independent behavior discards creep relaxation effects which enter eq. ([!]) in form of a 
diffusion term VfcVfc'I' with a rate independent coefficient. Such a term with an infinitesimally 
small prefactor is enough to regularizes the singularities and one gets the so-called viscosity 
solution, which we will consider in the following. 

We solve eq. (||) numerically on a finite difference grid using the problem solving environ- 
ment CACTUS 4.0 jll). The initial condition was a Gaussian random field with amplitude and 
correlation length one. The convolution was performed in Fourier space using the fast Fourier 
transform package fftw 2.1.3 We use an essentially non-oscillatory (ENO) scheme in 
combination with Godunov's method to minimize numerical damping p3| . The algorithm 
is proven to converge to the viscosity solution. The numerical grid has periodic boundary 
conditions; in one dimension (ID) the length of the system is 100 and the system had 4084 
points; in two dimensions (2D) the size of the system is 25 2 and the grid was 1016 2 . The 
three-dimensional (3D) simulation cell was 12. 5 3 and the grid size 264 3 . The ENO stencil 
width was 4 in one and 5 in 2D and 3D, respectively, and the time stepping scheme is a 
simple explicit Euler scheme with a time step 5a = 10~ 3 . The system is cycled by unloading 
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Fig. 2 - Blowup of the region 18 < x < 20 from our ID simulation. The solid lines represent 'i'(x) 
for a values 0.1 apart starting at a — 0, with curves further up having larger a. The dashed lines are 
'iz(x) while unloading and the dotted ones while reloading. The dotted lines are hardly visible because 
in the reversible region they lie on top of the dashed ones. At the cusps, the spatial discretization is 
visible. The inset shows the initial (solid) and final (dashed) configuration of the whole system, i.e. 
for a — and 6. 



beginning at various points a n and reloading at somewhat above zero stress (to reduce the 
effect of numerical damping) . 

Fig. ^ shows the evolution of the Burgers equation in ID, corresponding to eq. (||) with 
Sij — 1. The cusps develop in the valleys at a ~ 0.3, and the unloading and reloading paths 
are both shown. Since the cusps disappear immediately upon unloading, the order parameter 
evolution is reversible until the stress grows to match the previous maximum, so our model 
exhibits work hardening. 

Finally, Fig. |^ shows the cusp locations in two cross-sections of a 3D simulation with 
S X x = 2/3, S yy — S zz = —1/3, and zero otherwise, appropriate to a tension test. We 
see that the cusps form nearly flat ID interfaces separating cell-like volumes. While the cut 
parallel to the loading direction (a) shows the expected asymmetry, the cut perpendicular 
to the loading direction (b) shows an isotropic hexagonal cell structure. Simulations under 
shear in 2D (S xx = 1/2, S yy — —1/2, and zero otherwise) show similar cusps to Fig. 
These morphologies are reminiscent of cell structures formed in hardened fee metals 
particular the orientation of the cell walls is the same as observed in experiments, because the 
symmetry if the loading process is reflected in eq. ([!]). This, however, is only true for isotropic 
initial conditions. 
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Fig. 3 - Location of the singularities for tension test at a = 1 in (a) the xy-plsme and (b) yz-plane, 
that is parallel and perpendicular to the loading direction, respectively. The lines are contour lines 
SijViVjty = 1.5, 3, 6, 9, ... . Our symmetry analysis does not provide us a physical interpretation 
for 

The size distribution of the cells depends on the initial conditions and the cells coarsen 
in contrast to the cell refinement observed in experiments. Moreover, the cell walls in our 
model blur upon unloading, (see Fig. |^), much more than is seen experimentally. It is clear 
that these effects call for a more realistic extension of the present model, but we believe 
that the emergence of spatial structures, coupled with finite stress singularities which lead to 
irreversibility, is a generic feature which will be shared by a more elaborated non-linear field 
description of plasticity. 

We now wish to calculate a scalar, e.g., the total plastic strain, from the W field which 
illustrates the transition from reversible to irreversible behavior in a similar way as the stress- 
strain relation in Fig. (jl|). The term we are going to calculate and which we will call strain is 

^ffS M ((V fc ¥) (Vi*) V 2 *) . (3) 

The angle brackets are spacial averages. There are two reasons for focusing on this term. 

First it is sensitive to the formation of cusp singularities. At the cusp, the second derivative 
gives a delta function. Hence, the appearance of singularities will greatly affect the strain. 

Second, the term in eq. (^) is part of a gradient expansion. We assume (as common among 
engineers) that the deformation is in the direction of the applied stress Sij . If we expand the 
strain rate de/dt to second order in and fourth order in gradients the term in eq. (^) is 
one of the terms. Other terms are either insensitive to the singularities or they are surface 
terms (which would give zero average strain since we use periodic boundary conditions). 
There is, however, one extra term in such an expansion, which is highly singular, namely 
Sij (VfcVf'F) V 2 ^. This term would lead to a delta-function squared at a cusp, which we 
again neglect. 

Fig. |^ shows the resulting stress-strain curves for our plasticity theory. In ID and 2D the 
term in angle brackets in eq. fcfy is a total divergence: it is proportional to V^V^ 1 ]/) 3 and 
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Fig. 4 - Stress-strain relation according to eq. (g) for our ID (solid line) and 2D (dashed line) model, 
and for our 2D model with a stress-dependent prefactor 1 + 0.25 S mn S mn (inset). 



V a [(V x #) 3 /3 - (V^)(V y *) 2 ] - V y [(V^) 3 /3 - (V y f)(V^) 2 ] in ID and 2D, respectively. 
This implies that (dsij/dt) = 0, explaining why the unloading curves are vertical in Fig. [| 
Plastic deformation occurs on the main loading curve because the cusps act as sources for the 
total divergence term. (We do not show stress-strain curves for the 3D simulations. Because 
of numerical limitations the results would be governed by finite size effects, i.e., too few cells 
in the sample volume. Qualitatively the stress-strain curves look similar.) 

In the inset of Fig. || we added a stress dependent prefactor to cq. (||) in order to change 
the curvature of the stress-strain curve and make it look more similar to Fig. ([!]). This is 
to demonstrate that the upward-bending of the stress-strain curve is more related to the 
interpretation of the field \& rather than to the evolution equation (|l|) . The form of the initial 
conditions could also influence the shape of the stress strain curve considerably, in particular 
the onset of instabilities. 

Note that the point where the unloading curves meet the main stress-strain curve moves 
to higher stresses for increasing plastic deformation, i.e., the yield surface moves to higher 
stresses for higher deformation. In other words, our model equation shows work hardening. 
The work hardening is generic to the evolution equation and independent of the formula for 
the plastic strain, since the yield surface (i.e., the formation of singularities) moves to higher 
and higher stresses in eq. (|l|). 

The simple Hamilton- Jacobi equation (or anisotropic Burgers equation) (|l|) together with 
the plastic strain as calculated from eq. (^) forms a theory which generically has all three 
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mysterious features of plasticity, namely an evolving yield surfaces, cell structure formation 
(with the proper cell morphology for the given loading), and work hardening. However, as 
already mentioned, there are significant differences between our model and real plasticity. The 
cell structure in our model coarsens while it refines in real systems; also, the cell walls blur 
upon unloading in our model. Some of these shortcomings might be related to the choice of 
the viscosity solution as weak solution, to the absence of noise in our evolution equation, or 
the choice of initial conditions. 

The main point we want to make here is, that an evolving yield surface and work hardcning- 
like behavior is generated gencrically by the formation of cusp-singularities in Hamilton- Jacobi 
equations (and probably other types of singularities in different evolution equations). The 
analogy of the singularities in our model with the cell walls is only qualitative, but it suggests 
that a mesoscale continuum theory will have cell walls as singular structures, that will be 
central to work hardening and an evolving yield surface. We hope that the ideas presented in 
this paper will motivate further research in that direction. 
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